
%{
epsilonw      = 0.1;   % Learning rate for weights 
epsilonvb     = 0.1;   % Learning rate for biases of visible units 
epsilonhb     = 0.1;   % Learning rate for biases of hidden units 
weightcost  = 0.0002;   
initialmomentum  = 0.5;
finalmomentum    = 0.9;
%}


[numcases numdims numbatches]=size(batchdata);

if restart ==1
    restart=0;
    epoch=1;
    
    % Initializing symmetric weights and biases. 
    vishid     = 0.01*randn(numdims, numhid);
    hidbiases  = zeros(1,numhid);
    visbiases  = zeros(1,numdims);

    poshidprobs = zeros(numcases,numhid);
    neghidprobs = zeros(numcases,numhid);
    posprods    = zeros(numdims,numhid);
    negprods    = zeros(numdims,numhid);
    vishidinc  = zeros(numdims,numhid);
    hidbiasinc = zeros(1,numhid);
    visbiasinc = zeros(1,numdims);
end

error_plot = zeros(1,maxepoch);
tic
for epoch = epoch:maxepoch,
    errsum=0;
    
    rand_index = randperm(numbatches);
    
    for batch = 1:numbatches,
       
        %%%%%%%%% START POSITIVE PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        data = batchdata(:,:,rand_index(1,batch));
        flag_rate_template = data > 0;
        poshidprobs = 1./(1 + exp(-data*vishid - repmat(hidbiases,numcases,1)));    
        
        posprods    = data' * poshidprobs;
        poshidact   = sum(poshidprobs);
        posvisact = sum(data);
    
        %%%%%%%%% END OF POSITIVE PHASE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        poshidstates = poshidprobs > rand(numcases,numhid);

        %%%%%%%%% START NEGATIVE PHASE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        negdata = 1./(1 + exp(-poshidstates*vishid' - repmat(visbiases,numcases,1)));
        negdata = negdata.*flag_rate_template;
  
        neghidprobs = 1./(1 + exp(-negdata*vishid - repmat(hidbiases,numcases,1)));  
        
        negprods  = negdata'*neghidprobs;
        neghidact = sum(neghidprobs);
        negvisact = sum(negdata); 

        %%%%%%%%% END OF NEGATIVE PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        err= sum(sum( (data-negdata).^2 ));
        errsum = err + errsum;

        if epoch>5,
            momentum=finalmomentum;
        else
            momentum=initialmomentum;
        end;

        %%%%%%%%% UPDATE WEIGHTS AND BIASES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        if batch == numbatches
            count_last = sum(data')>0;
            vishidinc = momentum*vishidinc + ...
                epsilonw*( (posprods-negprods)/sum(count_last) - weightcost*vishid);
            visbiasinc = momentum*visbiasinc + (epsilonvb/sum(count_last))*(posvisact-negvisact);
            hidbiasinc = momentum*hidbiasinc + (epsilonhb/sum(count_last))*(poshidact-neghidact);
        else
            vishidinc = momentum*vishidinc + ...
                epsilonw*( (posprods-negprods)/numcases - weightcost*vishid);
            visbiasinc = momentum*visbiasinc + (epsilonvb/numcases)*(posvisact-negvisact);
            hidbiasinc = momentum*hidbiasinc + (epsilonhb/numcases)*(poshidact-neghidact);
        end
        vishid = vishid + vishidinc;
        visbiases = visbiases + visbiasinc;
        hidbiases = hidbiases + hidbiasinc;

        %%%%%%%%%%%%%%%% END OF UPDATES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    end
   % if mod(epoch,5)==0
        MAE
        toc
        tic
   % end
    %fprintf(1, 'epoch %4i error %6.1f  \n', epoch, errsum);
    error_plot(1,epoch) = errsum;
end;

plot(error_plot);